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Section  I 

PLASMA  TI1LORY 

There  are  no  additions  to  theory  to  report  this  quarter. 

Section  1 Z 
SIMULATION 

A.  DRIFT  CYCLOTRON  INSTABILITY  PARTICLE  SIMULATIONS 
Joe  Koo  Lee  (Prof.  C.  K.  Bird  sail) 

The  drift  cyclotron  instability  ( DC  I ) may  occur  in  any  finite 
co 1 I i s i on  I ess  magnetized  plasma  with  a Maxwellian  velocity  distribution. 

The  DCI  may  be  a major  confinement  problem  in  most  mirror  machines,  like 
2X1  IB  or  TMX , even  after  the  DCLC  mode  is  completely  stabilized. 

The  many  theoretical  studies  of  DCI  have  not  been  very  suc- 
cessful in  providing  uncontroversia I information  especially  on  the  nonlinear 


evolution  of  this  instability.  We  are  using  particle  simulations  to  obtain 
linear  and  nonlinear  behavior  of  this  instability. 

We  have  simulated  DCI  using  an  electrostatic  2 -d i mens i ona I parti- 

1 

cle  code,  E20HAR,  for  a Maxwellian  velocity  distribution,  f (v  1,  with  a 

density  gradient  in  a uniform  magnetic  field.  We  use  a slab  model  where 

the  plasma  density  changes  along  only  the  x direction,  and  the  magnetic 

field  is  along  the  z direction  as  in  Fig.  1.  The  temperature  is  uniform, 

~T  =0  (s  = e,i)  and  the  elect  rons  are  cold  initially  T =0  ; T.^0  . 
s e i 

The  densities  vary  as 


2 


Iwi th  A.  , A chosen  to  make  n =n.  . The  system  is  periodic  in  y , 

i e eo  i o 

bounded  in  the  outer  x , and  has  an  inversion  symmetry  at  the  left  x 
boundary  (the  zy  plane). 

As  reported  in  the  last  QPR,  our  simulations  produced  linear 

growth  rates  in  accordance  with  the  ROOTS  nonlocal  theory.  We  will  show 

detailed  results  for  two  mass  ratios,  namely,  m./m  =25  and  100  , in  an 

i e 

ERL  Report  in  the  next  quarter. 

4 The  nonlinear  evolution  of  this  instability  has  been  found  to 

be  complicated  in  single-mode  simulations.  We  observed  that  the  electro- 
static field  energy  saturates  at  1 'v 2 % of  the  initial  ion  kinetic 
energy  for  most  single-mode  cases  (of  about  20  runs).  The  possible  satur- 
ation mechanisms  may  be  listed  as: 


(1)  particle  (electron  or  ion)  trapping, 

(2)  nonlinear  shift  of  the  ion  cyclotron  frequency, 

(3)  wavebreaking  (of  electron  or  ion  wave), 

(A)  electron  (collisionless)  heating,  or 

(5)  the  density  gradient  modification  due  to  x B drift. 

Some  of  these  mechanisms  may  be  correlated  with  each  other.  We  will  present 


I 


B.  LOWER-HYBRID  DRIFT  INSTABILITY  SIMULATIONS  USING  AN  ESI  HYBRID  CODE 
Yu-Jiuan  Chen  (Prof.  C.  K.  Birdsall) 

In  the  last  QPR,  according  to  Eqs.  (3)  and  (A),  the  ions  can  be 

simply  treated  as  unmagnet  i .ted  particles.  Assume  ions  are  e lect ros tat i - 
ca  1 1 y conf ined,  i.e.,  v . = v_  + vv  = 0 where  v_  is  the  E x B drift  and 

-y,  - 1 - E - 

v y is  the  density  gradient  drift.  Therefore,  the  "ghost"  ion  method1 
mentioned  in  the  last  report  is  not  needed  to  simulate  the  lower-hybrid 
drift  instability  in  the  ions'  frame. 

The  electrons  are  treated  as  a linearized  fluid  with  suscep- 
tibility x (k,ui)  to  study  the  lower-hybrid  drift  instability.  This 
e 

method  was  originally  used  by  Bruce  Cohen  and  Gary  Smith.1  \^(k,u>)  is 
given  by 


.2  I - I (b)e  k k v,(k) 

xe  (k,w)  = + TIT  L — 

u)  b kvu)-kv,_ 

ce  y e y E 


for  a slab  shape*',  where 


b = k 


2 

' y m u» 
e ce 


v . (k  ) = 
" Y 


I (b ) e 

m u.  o 
e ce 


-b  vn  1 


v^  is  tne  E * B drift  velocity,  and 


and  T are  the 


•Bruce  I.  Cohen  and  Gary  R.  Smith,  "Efficient  Method  for  Local  Simulation  of 
Drift-Wave  Instabilities",  Proceedings  of  t‘*e  Ei_~.:~  Con*.  c-“  Numerical 
Simulation  of  Plasma,  PD -8,  1 S 7 8 . 

-R.  C.  Davidson  and  N.  T.  Gladd,  "Anomalous  Transport  Properties  Associated 
with  the  Lower-Hvbr id-Dr i ft  Instability",  Phvsics  of  Fluids.  Vol . 1 6 , No.  10, 
P.  1327  (1575). 
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electron  plasma  frequency,  cyclotron  frequency,  mass  and  temperature,  res 
pect i ve 1 y . 

The  Id  Fourier  transformed  Poisson  equation, 


[1  + xfi  (k,co)  ]<f>k  ( t ) = ~2  Pk(t) 


(3) 


may  be  written  as 


1 + 


pe 


ce/ 


(b) 


- b 


kvEW  + 

' k v 

e 


V-  /.  a 

k 


. ( i t- — kv_  ) p. 
2 \ at  E / k 


(4) 


by  replacing  w with  the  operator  i — , where  and  are  the 

perturbed  potential  and  ion  charge  density.  Put 


A = 1 


and 


W 1 

. -f  - 


(b). 


-b 


ce 


(5) 


8 £ 


kv  p 

pe  , . E 

^ KV  " A 


O 2 2 

2k  v 

e 


(6) 


Ea.  (4;  is  reduced  to  the  form 


T T *26;k  * "7  ( 1 nf  ‘ kV‘k)' 


(7) 


% 


— ■ - ■ 


_ 


_ 


Using  finite  differences,  we  obtain 


r 


4 


A 


0 - 


kvr 


A + B' 


C5  = KSQI 


(16) 


A + B' 


and 


kvi 


A ——  - B ' 

C6  E KSQI  • -TT"-2— 


(17) 


A + B1 


Finally,  we  obtain 


C1  = (ci  + ic2hk+  (c3+  iCA)pk+1  + (‘C5+iC6)pk  ’ 


(18) 


or 


Re  *"+1  = C1  Re  4>£  - C2  Im  + Z^  Re  p£+1  - Im  p"+1 


- Cr  Re  p"  - C,  Im  p" 
5 k 6 k 


(19) 


,n+  1 


Im  i.  = C 


Im  + C.  Re 
k 2 


*"  + C 


Im  P?+1  + C,  Re  P, 
3 k i*  k 


n+1 


- Cc  Im  pP  + C,  Re  . 
5 k 6 k 


(20) 


This  code  has  been  tried  in  the  1 ow-d r i f t- ve 1 oc  i t y regime 


2 2 2 2 

(v_<v.).  pa  ramete  rs  used  are  - /*  = 1 , - /-  .=n./m  =900 

Ei  pe  ce  pe  d i i e 


- . .:.t  = o.ok 

p I 


N.  = 102k  , v . = 0 . 1 A 1 A2  1 , v =0.QA  , v =0.u 

ion  i 


T ne 


simulation  result  is  shown  in  Fig.  1 


E ---  ■ e 
The  growth  rate  is  on  the  order 


— 


J 


r 


- 8 - 


of  the  theoretical  value  of  Davdison^.  B.  1.  Cohen  (LLL)  suggests  use 

2 

of  smaller  value  of  (u>  ./u)  .)  , and  a larger  number  of  particles,  due 

pi  ci 

to  problems  with  multibeaming  (numerical)  instabilities;  we  will  try  these. 


j] 

1 


9 


< 
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C.  FIELD  REVERSED  LAYER  SIMULATIONS  IN  2d 

Douglas  Harried  (Prof.  C.  K.  Birdsall,  Alex  Friedman) 

During  the  past  quarter  we  have  been  studying  the  two-dimen- 
sional electrostatic  particle  code  EZOHAR  with  a view  toward  modifying  it 
to  simulate  field-reversed  layers  (e.g.,  field-reversed  mirror  or  field- 
reversing  ion  layer  configurations).  The  simulations  would  thus  be  two- 
dimensional  across  the  r-theta  plane  of  the  cylindrical  field-reversed 
geometry.  This  type  of  simulation  should  enable  us  to  examine  the  theory 
of  Lovelace  *,  as  well  as  that  of  Uhm  and  Davidson2,  concerning  the  stab- 
ility of  long  layers.  When  electron  polarization  effects  are  included 
(which  will  not  be  in  the  early  stages  of  modification)  it  should  also  be 
possible  to  examine  Ogden's3  theory  concerning  the  microstability  of  these 
systems,  particularly  the  lower  hybrid  drift  instability. 

The  theories  of  interest  here  are  concerned  with  low  frequency  modes 
such  as  precession,  for  which  low-frequency  e lec t romagne t i c effects  are 
important.  Therefore  our  intent  is  to  proceed  by  converting  EZOHAR  from  an 
e I ec t ros t a t i c to  a Darwin  model  (i.e.,  neglecting  radiation  fields). 
Additionally  the  present  particle  code  will  be  converted  to  a quasineutral 
hybrid  one  in  which  ions  will  be  treated  as  particles  and  electrons  will 
be  treated  as  a fluid  (in  the  initial  modification,  a massless  fluid).  This 
should  be  adequate  for  the  frequencies  of  interest.  The  code  will  be 


j ‘R.  V.  Lovelace,  "Precession  and  Kink  Motion  of  Long  Astron  Lavers',  sub- 

mitted to  Physics  of  Fluids  (1978). 

‘H.  S.  Uhm,  R.  C.  Davidson,  "Stability  Properties  oF  a F i e 1 d- Reve rsed  Ion 
Laver  in  a Background  plasma",  PFC/JA-78-7  (1978''. 

■J.  M.  Ogden,  "Analytic  and  Numerical  St-dies  of  Microinstabi 1 i t ies  Driven 
by  Ions".  Ph.D.  dissertation,  Univ.  of  Maryland  (19781. 

i 


nonlinear,  and  serve  to  complement  linearized  codes  currently  being  developed 
(e.g.,  Friedman's4  three  dimensional  hybrid  code).  Allowing  nonlinearity 
may  complicate  observing  and  measuring  linear  growth  rates  between  noisy 
starts  and  saturation;  however,  problems  due  to  coarseness  of  initial 
conditions  usually  saturate  at  a low  level  in  a nonlinear  code,  rather 
than  grow  indefinitely  as  in  a linearized  system.  In  addition,  it  should 
be  of  interest  to  observe  saturation  and  other  nonlinear  effects  on  flute 
modes.  This  type  of  simulation  is  intended  to  provide  valuable  insight 
concerning  the  stability  of  f i e 1 d- reversed  layers. 


"A.  Friedman,  R.  N.  Sudan,  J.  Denavit,  "A  Linearized  Hybrid  Code  for  Stability 
Studies  of  Ax i svmriet r i c F i e 1 d- Re ve rsed  Equilibria",  Proceedings  of  the  Eightn 
Conference  on  the  Numerical  Simulation  of  Plasmas,  uor>tere\  CA,  June  1976- 
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Section  III 

CODE  DEVELOPMENT  and  MAINTENANCE 


ESI  ON  THE  CRAY 
Stephen  Au-yeung 


Dial  up  number: 


Same  as  the  A-ntachine 


Logging  in: 


c[user  no]  a 8 1 Os b 1 (ctrl-z) [combo] 


Getting  ESI: 


rf i ! em  / t v 

-read  1222  .cray  esI/C 


Extracting  files  from  esl/t: 

1 i b es l/£  / t v 

ok.  x esl/s  esl  esl/b/r 

ok.  end 


whe  re 


esl/s  is  the  FORTRAN  source 

esl  is  the  executable  file  and 

esl/b/T  is  the  relocatable  file  library 

Mod i fy i ng  ESI: 

qedesl/s  / t v 

- (qed  commands) 

- end 

Compiling  and  loading  ESI: 

eft  i = esl/s, b = esl/b,  = C /tv 

ldr  i = es1/b,x  = esl, lib  = es 1 /b/  - : t vS" 1 i b 
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7 . I nput  file  for  ES 1 : 

The  name  of  the  el  input  file  is  the  string  "input" 
followed  by  the  current  suffix  (channel).  For  example,  if 
suffix  b is  currently  used,  the  name  of  the  input  file 
must  be  "inputb".  This  allows  the  capability  of  running 
up  to  5 es 1 programs  simultaneously. 

The  input  file  specifies  the  problem,  plots  desired, 
etc.  For  example,  a two-stream  instability  problem  input 
could  be: 

box  b22  ee  stream  electron-electron  two  stream 
nsp  = 2 i rho  = 50  iphi  = 50  i xvx  = 25 
mplot  = 1 2 3 nt  = 300 
$ 

vO  = 1 . mode  =1  xl  = .001 
$ 

vO  = -1  . 

$ 

For  the  meaning  of  those  variables,  see  the  INPUT  section 
of  the  ESI  document  written  by  Bruce  Langdon. 

8.  Getting  hardcopy  output: 

In  order  to  get  graphics  from  the  nonimpact  printer 
one  has  to  first  transfer  his  files  to  the  A-machine.  The 
procedure  is  as  follows: 

netout  a alwith.  F 1 05  b. 

(crtl-d)  — to  logoff  the  Crav 

(login  under  the  A-machine) 

nipit  cfam.  filename  box  bnn  id  /tv 

Note  that  no  soace  in  "id"  is  allowed  and  if  the  error  mes- 
sage "The  F0  ...  file  is  missing"  appears,  ignore  it. 
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B.  EMI  Code 

No  special  report. 

C.  EZOHAR  DEVELOPMENT 

W.  M.  Nevins  and  Y.  Matsuda 

We  have  tested  EZOHAR,  a bounded,  2d3v  electrostatic  plasma 
simulation  code,  by  simulating  a plasma  slab  in  a uniform  magnetic  field. 
This  slab  is  uniform  in  y . The  plasma  density  is  constant  between  -L 
and  , and  zero  for  x > • This  configuration  was  chosen  for  test- 

ing the  code  because  the  linear  theory  is  relatively  simple.  Our  goal 
in  these  tests  are: 

(1)  To  compare  the  normal  mode  frequencies  with  the  linear  theory, 
as  this  provides  a test  of  the  accuracy  of  the  code  in  re- 
producing physics; 

(2)  To  test  the  power  of  EZOHAR  diagnostic  package.  In  par- 
ticular, we  wish  to  demonstrate  that  the  spatial  structure 
of  the  normal  modes  can  be  obtained  by  examining  the  fluc- 
tuation spectra  of  a stable  plasma. 

We  will  first  review  the  linear  theory  of  the  bounded  plasma 
slab,  and  then  present  a summary  of  our  simulation  results  to  date. 

A.  THE  MOMAL  ‘WES  OF  A BCUNVEV  PLASMA  SLAB 

Within  the  plasma  the  normal  modes  varv  spatially  like  sines 

and  cosines.  Hence,  there  are  four  linearly  independent  modes  for  each 

choice  of  k and  k . we  find  it  convenient  to  choose  these  four  modes 
x y 

to  be  eigenfunctions  of  the  inversion  operator,  and  write  the~  as 


i 

i 
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$j(x,y,t)  = 4>Jn  cos  k^x  cos  k^y  cos  (u»t  + <J> ^ ! 


:2(x,y,t)  = sin  k^x  sin  k^y  cos  (a>t  + <^2 ) 


?j(x,y,t)  = ^ sin  k^x  cos  k^y  cos  (u)t  + 4^) 


4^(x,y,t)  = i^n  cos  k^x  sin  k y cos  (iot  + 4^)  . 


Of  these  four  wave  forms,  only  and  $2  are  even  under  inversion,  i.e. 

4, (-x,-y,t)  = 4 1 (x.y.t)  (5) 


4>2(-x,-y,t)  - 1 2 ( x , y , t ) 


Hence,  the  inversion  symmetry  boundary  condition  imposed  at  x=0  (see 
QPR  of  July  1,  1976)  selects  out  modes  ij  and  4 , and  rules  out  modes 

4^  and  4^  , which  are  odd  under  inversion. 

For  x > , the  normal  modes  must  decay  exponentially  in  x 


4 . (x,  y , t ) = 4 U exp  (-k  x)  cos  k y cos  Lt  + :,) 

v v 


4,  (x,y,t)  = 


Out  / 

4-  exp  (-k  xj  sin  k y cos  Ut  + :,) 
2 y y 2 


We  w i 1 1 call  i ^ the  cos i ne- 1 i ke  mode , and  - the  "sine- 
like mode.  Matching  ; and  ?J/Vx  at  the  boundaries,  we  find  that 


r 


/ 

\ 


« 
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the  cosine- like  modes  obey 

tan  k L = k /k  (9) 

xx  y x 

while  the  sine-like  mode  obey 

cot  k L = -k  /k  . (10) 

xx  y x 

Periodic  boundary  conditions  are  imposed  in  y . Hence, 

k satisfies 

y 

k=^-m  , m = 0 , 1 , 2 . . . (ll) 

y y 

Specializing  to  the  case  , we  find  that  cosine- like 

modes  satisfy 

z tan  z = 2irm  (12) 


while  the  sine-like  modes  satisfy 

z cot  z = - 2irm  , m=  1,2,3  •••  (13) 


where 

a = k L . (14) 

x x 

Hence,  the  allowed  values  of  k may  be  labeled  by  two 
integers:  m is  the  number  of  nodes  across  the  system  in  y , and  n 

is  the  number  of  nodes  across  the  system  in  x . 

With  the  assistance  of  Miss  Jean  Yao  we  have  solved  Eqs. 
(12)  and  (13)  numerically,  with  a rootsolver  that  emplovs  the  global 
*'  wton's  method  (see  QPR  of  Oct.  1,  1977).  These  roots,  together  with 


1 


« 
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Eqs.  (1)  and  (2),  provide  the  allowed  values  of  |k|  , which  are  shown 
in  Tables  1 and  II,  and  the  normal  mode  structure,  an  example  of  which 
is  shown  if  Fig.  1.  The  normal  mode  frequencies  may  then  be  obtained 
by  solving  the  uniform  plasma  dispersion  relation  at  these  allowed  values 
of  | k | . This  is  accomplished  with  Gerver's  program  UNIFORM1,  available 
from  LIBRIS  on  the  MFECC. 


B.  ANALYSIS  OF  SIMULATION  RESULTS;  THE  FREQUENCY  AND  STRUCTURE  OF 

THE  NORMAL  MOVES 

In  our  computer  simulation  runs  we  save  time  histories  of 
the  potential  at  fixed  values  of  x k^  . The  E20HAR  pos t -processor , 

ZED  (see  the  previous  QPR  on  ZED  as  used  in  ESI),  then  calculates  the 
spectral  density  as  a function  of  the  discrete  variables  k and  x , 

y 

as  well  as  the  (essentially  continuous)  variable  to  from  these  time 

histories.  We  identify  the  peaks  in  plots  of  the  spectral  density  vs 

to  as  the  normal  mode  frequencies,  while  the  amplitude  of  each  peak  is 

2 

proportional  to  j $(x,k  ,u) j . Hence,  we  may  determine  the  spatial 

variation  in  the  amplitude  of  the  normal  modes  by  plotting  the  spectral 

density  vs  x for  fixed  k^  and  to  . 

Figures  2a  and  2b  of  the  April  1978  QPR  show  typical  plots  of 

the  spectral  density  vs  to  . The  upper-hybrid  frequency  and  the  electron 

gyroharmon ics  are  clearly  visible. 

In  general,  we  have  found  that  the  fluctuation  is  dominated  by 

the  long  wavelength  modes.  Hence,  those  modes  labelled  with  small  values 

of  m and  n have  oreater  amplitudes  than  those  modes  with  larger 
•M.  J.  Gerver,  "ROOTS,  A Disoerslcn  Equation  Solver",  ERL  Memo  No.  UCB/ERL 

M77/27,  31  Oct.  1976. 


m and  n values.  We  have,  as  yet,  only  been  able  to  identify  n = 0 
mode  structures  by  examining  the  fluctuation  spectrum  of  a plasma  slab 
in  thermal  equilibrium. 

The  technique  described  here  may  be  easily  generalized  to  obtain 
the  spatial  variation  in  the  phase,  as  well  as  the  amplitude,  of  the 
normal  modes.  This  phase  information  is  obtained  from  the  cross-spectrum 
of  the  potential  at  different  values  of  x . It  should  be  easier  to  measure 
the  spatial  structure  of  an  unstable  mode  (which  will  grow  to  large  ampli- 
tude) than  it  is  to  pick  the  mode  structure  out  of  the  thermal  fluctua- 
tion spectrum.  Hence,  we  expect  that  this  technique  will  be  very  useful 
for  comparing  the  normal  mode  structure  observed  in  computer  simulations 
with  the  results  of  linear  normal  mode  ca 1 cu 1 a t i ons . 
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D.  ESI  + EFL  Code 

Jae  Koo  Lee  (Prof.  C.  K.  Birdsall) 

A journal  article  is  in  preparation  for  the  Joun.naZ  o{j 
Compula  tic  not  Pkiji  i ci> . 


E.  RJET  DEVELOPMENT 

Stephen  Au-Yeung  (Alex  Friedman) 

The  PDP-11 /DMA  Controller  — model  121,  the  driver  of  the 
Versatec  Printer,  arrived  in  early  November,  1978.  However,  while  trying 
to  install  the  controller,  we  found  that  there  is  no  DD-11  DK  Expanded 
Board  left  in  the  PDP-11/04  computer.  Furthermore,  the  RJET  software 
development  group  discovered  that  16K  words  of  additional  memory  was  needed 
for  the  controller.  Orders  for  the  above-mentioned  items  have  been  placed 
and  we  are  waiting  for  DEC  people  to  deliver  them. 

F.  RUNNING  THE  SAME  PROGRAM  SIMULTANEOUSLY  UNDER  DIFFERENT  CHANNELS 
Stephen  Au-Yeung 

Inserting  the  following  code  at  the  beginning  of  a program 
is  an  example  of  how  to  make  the  program  run  under  different  channels  simul- 
taneous ly : 

parameter  (nf  = 3) 
integer  fnames(nf) 

data  fnames/8r  +prog,  84  input,  8r  output/ 

call  setfile  (fnames,  nf) 

call  dropfile  (fnames(l)) 

call  open  (5,  fnames (2),  0,  len) 

call  create  (6,  fnames(3),  1,  -1) 


where  setfile  is  a subroutine  that  appends  the  current  suffix  to  all 
names  given  in  the  array  "fnames".  To  get  this  subroutine,  type: 

filem  rds  1222  .lib  flib  end  / tv 

for  the  a-machine  users,  or  type: 

rfilem  read  1222  .Cray  util/ b/ 1 end  /tv 

for  the  cray-1  users. 
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G.  ROOTS  IMPROVEMENT 

Bruce  Cohen  and  Neil  Maron,  LLL  (A.  Friedman  and  S.  j-yeung) 


A revised  version  of  the  HOOTS  program  has  been 
developed  by  Bruce  Cohen  and  Neil  Maron  at  LLL  with 
the  capability  of  running  a family  of  cases  without  the  necessity 
of  the  user  restarting  the  program  for  each  case.  To  access  the  new 
version,  t ype : 

f i 1 em  rds  3054  root  rollOa  xroots  cosout 

The  first  of  these  files  is  the  FORTRAN  source,  the  second  is  the 
executable  file,  and  the  last  is  a sample  cosmos  deck  for  output 
process i ng . 


To  execute  the  new  version,  type:  roots  / t v 

A sample  terminal  interaction  session  is  as  follows: 


CSFR 

x roo  Is  /'  S 1 . 5 

PROG 

type  box  and  number  (box  Inn): 

CSFR 

box  u22 

PROG 

i npu t : 

CSFR 

more=l  j ay=0  0 0 nspec=3  charg=l  . 1 . -1  . 

omp s q=9 . e p 

s — . 3 1 .31 

.31 

PROG 

continue  namelist  input 

CSFR 

amass=l  1 . .00054  137  dense=l  . 125  -. 125 

1 r 1 a rm= 1 . 

33333  0 . 

end 

PROG 

i npu t : 

CSFR 

more=0  eps=.29  .29  .29  end 

PROG 

a 1 1 clone 

The  user  sets  more=l  if  one  or  more  case 

s ) will  f o 1 

ow  t he 

current  case,  and  more=0  if  the  current  case  is  the  last.  Namelist 
input  i.,  "erminated  with  "end". 


The  file  cosout  is  a cosmos  deck  to  send  output  to  the  versatec 


and  fiche 

least  the 

dev 

box 

ice,  and  should 

number  must  be 

be  modified  1 
changed.  To 

1 o suit  t he  user 
proce  s s oil  t pu  t , 

s needs  — at 
type  : 

cosmos 

i=cosout  / t v 
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H.  RINGHYBRID  CODE;  LOW  DENSITY  PLASMA  CONSIDERATIONS 

Alex  Friedman 

A.  INTRODUCTION 

The  linearized,  3d  hybrid  code  RINGHYBRID1  was  originally 
developed  at  Cornell  University  for  the  study  of  f i e 1 d - reversed  ion 
ring  stability  in  a dense  cold  background  plasma.  It  is  currently  being 
supported  by  this  author  at  U.C.  Berkeley.  The  code  follows  the  develop- 
ment in  time  of  perturbations,  with  azimuthal  mode  number  C , about  an 
axisymmetric  zero-order  state.  The  equilibrium  current  is  purely  azi- 
muthal and  is  carried  entirely  by  the  energetic  ion  component,  which  is 
modeled  by  discrete  particles  in  the  code.  [To  zero  order  these  are  axisym- 
metric rings,  while  to  first  order  they  are  perturbed  in  shape  by  a dis- 
placement eexp  (\lQ)  .]  Alternately,  discrete  particles  which  are  sta- 
tionary to  zero  order  can  be  initialized  uniformly  in  space  and  the  fluid 
ion  component  dispensed  with.  In  addition,  the  cold  background  of  ions 
and  a complement  of  inertialess  cold  electrons,  of  density  such  that  equil- 

1 

ibrium  charge  neutrality  obtains,  are  described  by  fluid  equations. 

One  potential  application  of  the  code  is  to  the  study  of  field- 
reversed  mirror  plasmas.  However,  in  this  application  the  assumption  of 
the  existence  of  a uniform  dense,  cold  plasma  component  is  not  appropriate. 
Thus,  we  need  to  examine  the  behavior  of  the  RINGHYBRID  algorithm  when 
applied  to  equilibria  containing  regions  where  the  density  is  low,  and  to 
consider  modifications  of  the  algorithm.  We  first  discuss  the  present 

*A.  Friedman,  R.  N.  Sudan,  and  J.  Denavit,  Pyoc.  the  Eighth  Confae-Xence  on 

Numerical  Simulation  o { Pla&mai,  Article  PC-13,  June  28-30,  1978,  Monterey, 
Cal i forn i a. 
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algorithm,  and  the  nature  of  the  limitation  this  algorithm  imposes  on  the 
size  of  the  t'mestep.  A revision  to  the  algorithm  which  eliminates  this 
limitation  is  proposed.  Finally,  some  considerations  which  arise  when 
trying  to  generalize  the  mode  I to  one  allowing  regions  of  vacuum  are  pre- 
sented . 

8.  THE  RINGHVBRW  ALGORITHM;  CONVERGENCE  CRITERION  FOR  THE  ITERATIVE 
SOLUTION  OF  OHM'  i L Alt' 

In  this  section  we  describe  the  algorithm  presently  employed 

in  the  solution  of  Ohm's  Law,  and  the  restriction  imposed  by  this  algorithm 

on  the  allowable  timestep  size,  At  . Defining  dimensionless  code  vari- 

2 2 2 

ables  as  B°  = u)°.At  , n.  = uj  .Ax  /c  , with  similar  units  for  other  quan- 
c i i p i 

tities  (to  eliminate  At  , Ax  in  the  particle  pusher  and  to  simplify  the 
algebra),  omitting  superscripts  on  first-order  quantities,  and  letting  a 
tilde  ( ~ ) denote  a quantity  defined  at  the  previous  time  level,  Ohm's  Law 
becomes 


n.  E = ( B°/2 ) xVxVxE  + W , 


(D 


whe  re 


W*  ; (B°/2 ) x v x V x | + b°xj.  - B°x  V x 6 . 


(2) 


Here  J.  is  due  to  the  discrete  ion  component  (we  assume  no  fluid  ion 
component  present),  and  the  last  term  in  W is  obtained  from  Jdt  VxVxE 
The  present  version  of  RJNGHYBRID  solves  Eq-  (0  by  direct 
iteration;  dividing  by  n.  , E (in  the  left  hand  side)  is  set  equal  to 
the  right  hand  side,  which  is  itself  a weak  function  of  E . 


the  right  hand  side,  which  is  itself  a weak  function  of 


This  i terat ion 


20 


2 2 2 

converges  when  (ic°jAt/2)(k  c / to  . ) < 1 (approximately).  Assuming 

k Ax  = 7i  , in  the  dimensionless  units  employed  this  becomes 
max  ' 


n B°/2n.  < 1 . 

i 


(3) 


Note  that  the  Courant  condition  on  Alfven  waves  restricts  B°//n7  to 

i 

small  values,  so  that  as  n.  decreases  the  convergence  criterion  for  the 
direct  iteration  eventually  imposes  a more  severe  limit  on  At  than  does 
the  Courant  condition.  Inhomogeneous  mirror  plasmas  will  contain  regions 
of  small  n.  , so  that  an  undesirably  small  At  must  be  chosen  if  the  pre- 
sent algorithm  is  used. 

C.  A PROPOSED  MODIFICATION  TO  FACILITATE  INVERSION  OF  OHM' i LAW  IN 
PROBLEMS  INVOLVING  SMALL  n- 

t 

To  eliminate  the  restriction  imposed  by  direct  iteration  on  (1) 
it  is  necessary  to  bring  the  term  (B°/2 ) *V*V*E  to  the  left  hand  side. 
Since  B°  and  n.  are  functions  of  the  spatial  variable  z a direct  in- 

i 

version  by  Fourier  analysis  in  z and  tridiagonal  solution  in  r is 
not  feasible;  a form  of  ICCG  inversion  would  be  suitable,  but  for  sim- 
plicity a simpler  iterative  method  is  proposed.  This  would  replace  the 
direct  iteration  on  (l). 

Equation  (1)  represents  six  equations  for  the  unknown  components  of 

E , since  E , E , E each  have  2 parts,  EC  ® cos  £6  and  ES<*sin  £0  , 

~ r 0 z 

where  £ is  the  azimuthal  mode  number.  Thus,  in  each  cell  (i,j)  we  have 

a 6 >6  matrix  to  invert  for  E..  , assuming  E.  . . , E..,  . , ...  known. 

-ij  i,J+1  i+1,J 

For  example,  the  first  equation  (for  the  C component)  is,  letting  a = Ar/Az 


% 


n . E . . ■ \U/r.)C..  - 2 + 2cx+ l/r^Et.  . 

i r i j 2 l J r.j  j e i j J 


B° . . 


= W 


nj 


|(af/2rj 


)(E  -E  ) 

i+1,j  i-1,j 


- J(EC,  + ) - (1  + 1/2 r . ) E*T 


"0.  . . e.^.  . 

1-1, J 1 + 1, j 


j ei,j+1 


- (1  - 1 / 2 r . ) E^  + U/2r.)  (ES 

i,j-i 


- E' 


J ri,j+1 


XC.  . . 
«*i  j 


(4) 


Similarly,  the  other  5 component  equations  are  found,  yielding  an  equa- 
tion of  the  form 


The  "E.  . ”,  etc.  parts  of  X already  appear  in  the  current  version  of  RING 

- i , j + 1 


HYBRID,  and  only  the  "E..”  terms  need  be  deleted  from  the  source  terms  in 

x -i  j 


the  present  code.  Explicitly,  we  find 
letting  b * B? ./2 


A to  be,  omitting  indices  i ,j  and 


- 29 


This  matrix  is  singular  for  n.  =0  , and  ill-conditioned  when 

n.  is  very  small;  to  see  this  note  the  linear  dependence  of  the  first 

. c s 

and  fifth  rows  when  n.=0  . Inverting  this  matrix  gives  E ..  , E ..  , ... 

i r i j r i j 

as  a function  of  W . . and  the  EC which  are  assumed  "known" 

~ i J r i , j+1 , 

from  the  last  pass  in  the  iterative  scheme  proposed.  This  inversion  should 

— -1 

be  done  once,  on  paper,  and  A programmed  into  the  code. 

While  this  modification  does  not  change  the  model  equations  at 
all,  and  hence  in  principle  does  not  alter  the  range  of  problems  for  which 
the  model  equations  provide  a reasonable  description,  the  relaxation  of  the 
restriction  on  At  should  make  practical  the  simulation  of  inhomogeneous 


mi rror  p 1 asmas . 


V.  GENERALIZATION  TO  PROBLEMS  INVOLVING  VACUUM  REGIONS 


A possible  generalization  of  the  code  involves  matching  the 

plasma  model  (with  either  direct  iteration  or  the  above  method)  to  a vacuum 

2 

region,  where  the  equation  V E = 0 is  assumed  to  hold.  Note  that  speci- 


fication of  VxVxE  = 0 in  vacuum  is  by  itself  insufficient.  The  addi- 
tional spec  ificationof  V • E = 0 or  at  least  v(v  • E)  = 0 is  necessary  to 

determine  all  components  of  E . Combi  ning  V • E = 0 with  VxV*E  = 0 
2 

yields  V E = 0 as  the  equation  to  be  solved. 

The  author  has,  for  a previous  version  of  the  3d  linearized 

2 

code,  developed  an  iterative  solver  for  the  equation  V A = J ; with  a change 

2 

of  notation  this  is  directly  applicable  to  the  solution  of  V E = 0 . 
Explicitly,  in  "vacuum"  cells  (possibly  defined  to  be  those  cells  where, 
say,  v^j ^venAt/Ax > 0. 5) , we  solve 


S, . . - UCR . E° . . - UCX.E* . = 0 

1 ij  j ri j j 0i j 


S„.  . - UCR.E^. . - UCX.EC. . = 0 
2 i j j 0ij  j rij 


(7) 
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where  UCR^  = 2+2a2+£2/r2+1/r2  , UCX ^ ; 2£/r?  , and  come  from  cells 

of  the  unknown  other  than  i , j , and  from  a source  term  taken  to  be  zero 
for  the  current  application.  We  solve  to  find 


(UCR . S . . - UCX.S0.  ., 
J 1 ' J J 2i j 


(UCR? . - UCX2.; 
• J ' J 


Es.  . = UCR  UCX 
eij 


(8) 


and  similarly  for  the  other  components. 

However,  it  remains  to  be  determined  whether  the  solution  to 
this  equation  in  vacuum  implies  a smooth  transition  from  the  solution  to 
Eq . (l)  in  plasma.  In  a private  communication,  J.  Byers  suggests  that 
in  the  last  plasma  cells  or  first  vacuum  cells  strong  "sheath"  effects 
may  be  present.  Some  observations  relevant  to  this  issue  follow. 

Where  n.  = 0 , the  plasma  Eq . (1),  Ohm's  Law,  makes  little 

physical  sense.  However,  formally  it  reduces  to 


- B°  x v x V x (E  + E)/2  = - B°  x V x B , 


i .e.  , 


B°  x v x B = 0 , 


or  , 


B°  x v x v x E = 0 


(9) 


(10) 


(H) 


with  proper  choice  of  initial  conditions.  This  last  is  certainly  true,  and 
is  a "subset"  of  VxVxE=0  , but  evidently  is  not  sufficient  to  specify 
E in  vacuum  (even  when  n.  is  nonzero  but  very  small,  problems  are  likely 
to  be  present).  Consider  the  Id  , rad i a 1 -var i at i on-on ly  limit:  Eq.  (1) 


becomes 


! 

! 


4 

i 


r:  n . . E . + bz  . , 2 + 1 A ' 

•J  ri 


1 + \ 2r 


1-t  2 r . ) E ) (12a) 


n.  .E  = W . = B" 

‘J  "J  j 1 r 


(12b) 


The  laiterof  these  gives  E c\  . j ri . . , wni  cb  i s a bad  1 y posed  equa- 

tion when  n.  is  sma I 1 ince  J.  w i 1 : f1uctuate  due  to  the  small  number 
i i r 

of  part i c les/ce  I 1 (n.  i 5 assumed  "frozen",  as  it  is  a zero-order  quantity). 

With  E so  specified,  J ......  is  specified  and  hence  E is  deter- 

I u i A L r 

mined.  In  vacuum, 

(2+1/r2)E.  = (1  + 1/2r.)E...  + (l-l/2r.)E.  , , (13) 

J -J  J -J+1  J 'J*1 

the  same  equation  holding  for  both  components  E and  E . Even  though 

r o 

n.  may  be  small  in  (12a),  n.E  may  not  be,  and  in  any  event  E is 
i i r 0 

determined  by  ( 1 2 b ) . Thus,  the  vacuum  and  plasma  equations  do  not  appear 
at  all  similar;  only  the  latter  supports  waves,  and  the  interface  must 
reflect  these  waves  in  some  manner,  so  that  sheath  effects  are  likely. 
Appropriate  jump  conditions  at  the  interface  must  thr  be  formulated.  These 
questions  are  undergoing  further  study. 

E.  CONCLUSIONS 

In  its  present  form  the  RINGHYBRID  code  is  applicable  to  a 
restricted  class  of  problems  involving  dense  plasmas  which  are  spatially 
nearly  homogeneous.  A modification  has  been  proposed  which  would  facilitate 
modeling  of  more  strongly  inhomogeneous  mirror  plasmas,  as  in  f i e I d- reversed 
mirror  configurations.  Finally,  some  considerations  which  arise  when  one 
attempts  to  generalize  the  model  itself  to  include  vacuum  regions  have  been 


presented . 
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Section  IV 

PLASMA  SIMULATION  TLXT 

Our  800  pages  of  notes  have  been  edited  and  redone,  with  con- 
tents as  below.  The  new  chapters  6 , 7 , 1 2 , 1 3 , 1 4, 1 5 , 1 6 , and  the  Appendices. 

The  set  is  available  from  Prof.  Birdsall.  Send  a check  for  $23-75  (our 
cost)  made  out  to  Regents,  University  of  California.  The  printing  was 
limited;  first  come,  first  serve. 

PLASMA  PHYSICS  vta  COMPUTER  Sl-UIATION 
bj 

CI1J.1U4  K.  Siidsail  and  A.  S*uc*  lAjtgdon 


Tart  I PRIMER 

One  Dimensional  Electrostatic  and  Electromagnetic  Codes 

Chapter  1.  Some  ideas  on  why  computer  simulation  using  particles  rakes 
good  physical  serse. 

2.  Over-all  view  of  a Id  electrostatic  program. 

3.  A id  electrostatic  program,  £ Si  . , 

4 Introduction  to  the  numerical  methods  used. 

5.  Projects  for  ESI. 

6.  A Id  electromagnetic  program,  ENl . 

7 . Projects  for  EMI . 

Part  II  THEORY 

Plasma  Simulation  Using  Particles  in  Spatial  Grids 
with  Finite  Time  Steps;  Varr  Plasma 

8.  Effects  of  the  spatial  grid. 

9.  Effects  of  the  finite  time  step. 

10  Energy  conserving  simulation  models. 

11.  D1 pole  model s . 

12.  Kinetic  theory  for  fluctuationa  and  notae. 

13.  Statistical  mechanics  of  a sheet  plasma. 

Part  III  PRACTICE 

Programs  In  Two  and  Three  Dimensions;  Design  Procedures 

14.  Electrostatic  Programs  In  2d,  3d. 

15.  Electromagnetic  programs  In  2d,  3d. 

16.  Design  of  computer  experiments;  choice  of  parameters. 

Part  TV 
Appendices 

DECEMBER  1978 
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